clc;
clear all;
close all;
c=3.0e8;  
B=768e6;       %调频带宽
K=29.982e12;       %调频斜率
T=B/K;         %调频周期
Tc=160e-6;     %chirp总周期
fs=1e7;       %采样率
f0=77e9;       %初始频率
lambda=c/f0;   %雷达信号波长
d=lambda/2;    %天线阵列间距

N = 256;       %距离向FFT点数
M = 128;       %多普勒向FFT点数
Q = 180;       %角度FFT
%读取数据，能显示每帧数据
%%% This script is used to read the binary file produced by the DCA1000
%%% and Mmwave Studio
%%% Command to run in Matlab GUI -
%readDCA1000('<ADC capture bin file>') function [retVal] = readDCA1000(fileName)
%% global variables
% change based on sensor config
numADCSamples = 256; % number of ADC samples per chirp
numADCBits = 16; % number of ADC bits per sample
numRX = 4; % number of receivers
numLanes = 2; % do not change. number of lanes is always 2
isReal = 0; % set to 1 if real only data, 0 if complex data0
n_chirps=128;%每一帧中chirps数，no of chirp loops 
frame=50;%第多少帧
%% read file
% read .bin file
fileName='D:\20241022\4.bin';
fid = fopen(fileName,'r');
adcData = fread(fid, 'int16');
% if 12 or 14 bits ADC per sample compensate for sign extension
if numADCBits ~= 16
    l_max = 2^(numADCBits-1)-1;
    adcData(adcData > l_max) = adcData(adcData > l_max) - 2^numADCBits;
end
fclose(fid);
fileSize = size(adcData, 1);
% real data reshape, filesize = numADCSamples*numChirps
if isReal
    numChirps = fileSize/numADCSamples/numRX;
    LVDS = zeros(1, fileSize);
    %create column for each chirp
    LVDS = reshape(adcData, numADCSamples*numRX, numChirps);
    %each row is data from one chirp
    LVDS = LVDS.';
else
    % for complex data
    % filesize = 2 * numADCSamples*numChirps
    numChirps = fileSize/2/numADCSamples/numRX;%总chirps数
    numframe=numChirps/n_chirps;%总帧数 no of frames，总chirps数除以每帧chirps数。
    LVDS = zeros(1, fileSize/2);
    %combine real and imaginary part into complex data
    %read in file: 2I is followed by 2Q
    counter = 1;
    for i=1:4:fileSize-1
        LVDS(1,counter) = adcData(i) + sqrt(-1)*adcData(i+2); 
        LVDS(1,counter+1) = adcData(i+1)+sqrt(-1)*adcData(i+3); 
        counter = counter + 2;
    end
        % create column for each chirp
        LVDS = reshape(LVDS, numADCSamples*numRX, numChirps);
        %each row is data from one chirp
        LVDS = LVDS.';
end
%organize data per RX
adcData = zeros(numRX,numChirps*numADCSamples);
for row = 1:numRX
    for i = 1: numChirps
        adcData(row, (i-1)*numADCSamples+1:i*numADCSamples) = LVDS(i, (row-1)*numADCSamples+1:row*numADCSamples);
    end
end
% return receiver data
%adcData;接收的总数据retVal = adcData(:,(frame-1)*numADCSamples*n_chirps+1:frame*numADCSamples*n_chirps);%读取每帧数据
% 初始化framedata矩阵
numFrames = 100; % 总帧数
framedata = zeros(numRX, numADCSamples * n_chirps, numFrames); % 预分配三维矩阵

for frame = 1:numFrames
    framedata(:, :, frame) = getFrameData(frame, adcData, numADCSamples, n_chirps);
end
%framedata包含四个通道所有一百帧的数据
%下面是单独通道的数据
data_radar_1 = reshape(adcData(1,:),numADCSamples,numChirps);   %RX1
data_radar_2 = reshape(adcData(2,:),numADCSamples,numChirps);   %RX2
data_radar_3 = reshape(adcData(3,:),numADCSamples,numChirps);   %RX3
data_radar_4 = reshape(adcData(4,:),numADCSamples,numChirps);   %RX4
% data_radar=[];            
data_radar(:,:,1)=data_radar_1;     %三维雷达回波数据
data_radar(:,:,2)=data_radar_2;
data_radar(:,:,3)=data_radar_3;
data_radar(:,:,4)=data_radar_4;
%% 主程序

[range,speed,angle]  = threedfft(data_radar, N, M, Q, numADCSamples,numChirps,numRX);
figure;
speed1 = reshape(speed(:,:,1),N,M);   
speed_Temp = speed1';
[X,Y]=meshgrid((0:N-1)*fs*c/N/2/K,(-M/2:M/2-1)*lambda/Tc/M/2);
mesh(X,Y,(abs(speed_Temp))); 
xlabel('距离(m)');
ylabel('速度(m/s)');
zlabel('信号幅值');
title('2维FFT处理三维视图');
xlim([0 (N-1)*fs*c/N/2/K]);
ylim([(-M/2)*lambda/Tc/M/2 (M/2-1)*lambda/Tc/M/2]);

%% 函数
function [range,speed,angle]  = threedfft(data_radar, N, M, Q, numADCSamples,numChirps,numRX)
%距离FFT
range_win = hamming(numADCSamples);   %加海明窗
doppler_win = hamming(numChirps);
range_profile = [];
parfor k=1:numRX
   for m=1:numChirps
      temp=data_radar(:,m,k).*range_win;    %加窗函数
      temp_fft_1=fft(temp,N);    %对每个chirp做N点FFT
      range_profile(:,m,k)=temp_fft_1;
   end
end

%多普勒FFT
speed_profile = [];
parfor k=1:numRX
    for n=1:N
      temp=range_profile(n,:,k).*(doppler_win)';    
      temp_fft_2=fftshift(fft(temp,M));    %对rangeFFT结果进行M点FFT
      speed_profile(n,:,k)=temp_fft_2;  
    end
end

%角度FFT
angle_profile = [];
parfor n=1:N   
    for m=1:M   
      temp=speed_profile(n,m,:);    
      temp_fft_3=fftshift(fft(temp,Q));    %对2DFFT结果进行Q点FFT
      angle_profile(n,m,:)=temp_fft_3;  
    end
end
    %[range,speed,angle] = [range_profile,speed_profile,angle_profile]
    range = range_profile;
    speed = speed_profile;
    angle = angle_profile;
end

function retVal = getFrameData(frame, adcData, numADCSamples, n_chirps)
    % 检查输入参数
    if frame < 1
        error('帧数必须大于0');
    end
    
    % 计算帧的起始和结束索引
    startIdx = (frame - 1) * numADCSamples * n_chirps + 1;
    endIdx = frame * numADCSamples * n_chirps;
    
    % 检查索引是否超出数据范围
    if startIdx > size(adcData, 2) || endIdx > size(adcData, 2)
        error('帧数超出数据范围');
    end
    
    % 返回指定帧的数据
    retVal = adcData(:, startIdx:endIdx);
end

